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The pygmy dipole resonance in neutron-rich nuclei is studied within the framework of the Vlasov 
equation which is solved numerically. The interaction used in the Thomas-Fermi ground state and 
in the Vlasov equation is derived from an energy functional which correctly describes the equation 
of state of nuclear matter and neutron matter. It is found that the pygmy resonance appears in 
the electric dipole response of all nuclei with strong neutron excess, the energies and transition 
probabilities being in reasonable agreement with experimental results. Since the Vlasov equation 
does not account for any shell effects, this indicates that the existence of the pygmy resonance is 
a generic phenomenon and does not rely on the specific shell structure. Besides the electric dipole 
response, the isoscalar toroidal response is calculated. The transition densities and velocity fields are 
discussed. A comparison of the peak positions and velocity fields suggests that the pygmy resonance 
can be identified with one of the low-lying modes excited by the isoscalar toroidal operator. 



PACS numbers: 21. 10. Re, 24.30.Cz, 03.65. Sq, 02.70.Ns 

I. INTRODUCTION 

Neutron rich nuclei have become a very popular object 
of experimental and theoretical nuclear structure studies. 
Besides the crucial role these nuclei play in nuclear astro- 
physics and their importance for constraining the nuclear 
energy density functional, these nuclei exhibit fascinating 
properties which are qualitatively different from ordinary 
nuclei. For instance, as a consequence of the "neutron 
skin" surrounding the core of medium-mass and heavy 
neutron-rich nuclei, there are new kinds of collective mo- 
tion which are absent in nuclei without strong neutron 
excess (see [l[ for a recent review). 

A famous example for such a collective mode is the 
so-called pygmy resonance. Contrary to the well-known 
isovector giant-dipole resonance (GDR), where neutrons 
and protons move against each other, the pygmy-dipole 
resonance (PDR) consists, roughly speaking, of an oscil- 
lation of the neutron skin against the N ~ Z core. This 
mode is not only interesting in itself, but its existence has 
also a strong effect on the abundances of the elements 
in the universe Q. After first studies within schematic 
hydrodynamic models 0, H[ , the pygmy mode was inves- 
tigated within the random-phase approximation (RPA), 
using non-relativistic [f| or relativistic formalisms (T 
and beyond, using the quasiparticlc-phonon model [ 

Another exotic type of collective motion is the 
"toroidal dipole mode". This isoscalar mode, which is 
characterized by a velocity field of toroidal shape, was 
predicted many years ago . Since the restoring force 
for this kind of motion is generated by the distortion of 
the Fermi surface, this mode cannot be described within 
hydrodynamic models. Studies of this mode were car- 
ried out within the method of Wigner function moments 
[ill ], within nuclear fluid dynamics fl2| . and within the 
relativistic RPA [l3|]. Since the toroidal dipole mode is 
isoscalar and exists also in N — Z nuclei, it was usually 
not connected with the pygmy resonance, although dif- 
ferent calculations @, ||| snowed that the velocity field of 
the pygmy mode has a toroidal shape. 



In the past, semiclassical approaches such as the 
Steinwedel- Jensen model [l4| contributed a lot to the 
understanding of giant resonances and how they can be 
related to global properties of nuclei. Since semiclassi- 
cal approaches average over shell effects [l5j|, they lead 
to clear and intuitive pictures which, contrary to the re- 
sults of quantum mechanical RPA calculations, are not 
obscured by details of the specific single-particle ener- 
gies and wave functions of the nuclei under considera- 
tion. This is why new modes of excitation, such as the 
torus mode mentioned before, were often first identified 
in semiclassical approaches. 

In the case of the pygmy mode, a generally accepted 
picture is still missing. It is not even clear whether the 
pygmy mode is a generic collective mode like the giant 
resonances, or whether it depends on a particular struc- 
ture of the single-particle levels. The aim of the present 
paper is therefore to study this mode by solving the semi- 
classical Vlasov equation for the case of neutron-rich nu- 
clei. The Vlasov equation has been shown to give a rea- 
sonable description of generic properties of different col- 
lective modes of nuclei [l6l - fl8j . As we will see, the main 
result of the present study is that the pygmy mode is 
closely related to a low-lying isoscalar torus mode, which 
is clearly collective and whose existence is not limited to 
nuclei with neutron excess. 

In this paper, we are looking for a numerical solution of 
the full Vlasov equation without any additional simpli- 
fying assumptions, as opposed to, e.g., fluid-dynamical 
approaches. A solution of the Vlasov equation for col- 
lective modes was given in Refs. [lfl 0], but the calcu- 
lations were restricted to Woods-Saxon potentials with 
separable residual interactions. This is not sufficient for 
the description of modes which possibly depend on ex- 
otic ground state properties such as the neutron skin. In 
the present work, the starting point is the ground state 
in Thomas-Fermi (TF) approximation, calculated self- 
consistcntly with the same interaction that is also used 
in the Vlasov dynamics. The importance of a consis- 
tent description of the ground state and of the dynamics 
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within a transport model was recently pointed out in Rcf. 
[l9l | in the context of the monopole mode. 

For the calculation of the mean field entering both 
the TF and the Vlasov calculations, an effective interac- 
tion capable of describing exotic nuclei is needed. In the 
present work, a simplified version of the so-called BCP 
functional [201 ] will be used, which is an energy functional 
whose bulk part is based on a fit to microscopic Bruckner 
calculations, reproducing the equation of state of nuclear 
matter in a range of asymmetries from symmetric matter 
to pure neutron matter and in a range of densities from 
zero to more than saturation density. However, it seems 
unlikely that the general findings presented here depend 
on the details of the interaction. 

The paper is organized as follows. In Sec. m the 
method is explained. In Sec. IIII1 results for the elec- 
tric dipole response are presented. Besides the strength 
function, the transition densities and velocity fields of 
the GDR and the PDR are discussed. In Sec. IIV1 results 
for the response to the isoscalar toroidal dipole operator 
are shown. In particular, by comparing the transition 
densities and velocity fields with those discussed in Sec. 
IIII1 we see that in neutron-rich nuclei, the PDR and the 
GDR arc modes which are excited by both the electric 
dipole and the isoscalar toroidal dipole operators due to 
the mixing of isoscalar and isovector modes. Sec. [V] is 
devoted to the summary and conclusions. 



II. METHOD 

A. Vlasov equation 

Like other collective vibrations, the pygmy resonance 
has been studied within the random-phase approxi- 
mation (RPA) @-0], which can be interpreted as the 
small-amplitude limit of the time-dependent Hartree- 
Fock (TDHF) theory [ll|. Written in terms of the one- 
body density matrix p and the mean-field hamiltonian h, 
the TDHF equation reads 1 



iHp = 



■P 



(1) 



In the semiclassical h — > limit, the TDHF equation 
reduces to the Vlasov equation [ill [l8| ■ In order to see 
this, it is useful to work with the Wigner transforms of 
p and h, which are the distribution function /(r,p,t) 
and the classical mean- field hamiltonian h(r,p,t). For 
example, in the case of a purely local mean field U, the 
latter can be written as 

.2 



h(r,p,t) 



P 
2m 



+ U(r,t). 



(2) 



1 Protons and neutrons have of course different density matrices 
p p and p n , different mean-field hamiltonians h p and h„ , etc. In 
order to improve the readability, isospin indices a = p, n are 
omitted in this paper except when they cannot be avoided. 



Note that the mean field U depends on time through 
the time dependence of the density. To leading order 
in h, the Wigner transform of the commutator in Eq. 
(U} reduces to the Poisson bracket of the corresponding 
Wigner transforms, and one obtains the Vlasov equation 

f = {h f} = — ■ — - — ■— (3) 
' dv dp dp dr 

A discussion of the Vlasov equation as limiting case of 
more general transport equations can be found in Rcf. 



B. Numerical method 

In order to solve the Vlasov equation ([3]) numerically, 
we will employ the test-particle method which has of- 
ten been used for the description of heavy-ion collisions 
fj~8l . The basic idea of this method is to replace the 
distribution function f(r, p, t) by a finite number of delta 
functions ("test particles") 

MA 

f(r, P,t) = jfJ2 S ( F - r^MP - P* (*)) , (4) 

i—l 

where f\f denotes the number of test particles per nu- 
cleon and A is the mass number of the nucleus. In order 
to satisfy the Pauli principle, the density of test particles 
of each species (a = n,p) in phase-space must not ex- 
ceed 2Af / (2irh) 3 (the factor of 2 is the spin degeneracy). 
Inserting Eq. (|4]) into Eq. ([3]), one finds that each test 
particle has to follow its classical trajectory given by 



dp t 



(5) 



In the case of a purely local mean field, the equations of 
motion reduce to 



Pi 

rn 



VU{v u t). 



(0) 



It is clear that these equations of motion prevent the 
test particles from entering the classically fobidden re- 
gion. Note that this absence of tunneling is inherent to 
the Vlasov equation and independent of the numerical 
method. 

The density p(r, t) corresponding to the distribution 
function Q is a sum of delta functions and hence not 
suitable for any practical calculation. In order to obtain 
a well-defined density which can be used, e.g., for the cal- 
culation of the mean field U, it is common to replace the 
delta functions in Eq. (01 by Gaussians [22|, [23[ , leading 
to a smooth density 



PM)=]T 



AM (r _ r! ( t)) Vd 2 



(7) 



For the sake of consistency, if one uses the smooth 
density p(r) in the calculation of the mean field U(r), one 
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has to modify also the acceleration equation and replace 
the force at Yi by a force averaged over the Gaussian 

MM, 



where 



d 3 t 



'Id- 



U(r-s,t). 



(8) 



(9) 



Contrary to quantum molecular dynamics (QMD) [24[ 
and related approaches, where similar Gaussians (in r 
and p space) are used to simulate quantum effects, we 
wish to stay here in the semiclassical framework and 
therefore do not attach any profound meaning to the 
smoothed density p. We consider it as an auxiliary quan- 
tity one has to introduce in the test-particle approach in 
order to be able to calculate well-defined densities and 
mean fields. As we will see in the next subsection, an in- 
teresting aspect of the smoothing of densities and mean 
fields is that it acts exactly like a finite-range interaction. 



C. Interaction 

Until now, the mean field U entering the hamiltonian 
h has not been specified. On the one hand, it should 
be local for simplicity, and on the other hand, it should 
not be too simplistic if one wants to describe exotic nu- 
clei. In this work, it will be derived from the bulk part 
E^ t [p p , p n ] of the Barcelona-Catania-Paris (BCP) energy 
functional [20| . which is a parametrization of Bruckner 
G-matrix results for nuclear and neutron matter. If the 
interaction energy is written as 



E£k\Pp,Pn]= / d d r e£ t (p p , p n ) 



(10) 



the mean fields for protons (a = p) and neutrons (a = n) 
are given by 



a po 

U a ( Pp , Pn ) = -^ 

Op a 



(11) 



As explained in Sec. Ill B[ U is calculated from the 
smoothed densities p, i.e., 



U a (r,t) = U a {p p {r,t),p n (r,t)) 



(12) 



In addition to the bulk part E^ t , the BCP functional 
contains a finite-range part E™, a spin-orbit part E s '°', 
and a Coulomb part Eg- The finite-range part, which 
was introduced in Ref. [20[ in order to get the right sur- 
face energy, has mainly the effect to smooth out the mean 
fields as compared to the densities. The same effect can 
be achieved without any additional term if one intro- 
duces a finite range into the bulk term E-£ t . For instance, 
in the latest version of the BCP functional [H, has 
been substituted by a finite range in the quadratic term 
of ES°,. 



Here, instead of implementing directly a finite-range 
force between the test particles, we first calculate the 
smooth density by folding the distribution function with 
a Gaussian of width d [Eq. (J7J]. From this smoothed den- 
sity, we calculate the mean field [Eq. (fl2|) ] which is folded 
again by a Gaussian of width d when the force on a par- 
ticle is calculated [Eq. ©]. If ej^ t was quadratic in the 
densities, this procedure would be equivalent to a Gaus- 
sian finite-range interaction with width y/2d. Hence, if 
we choose d = 0.7 fm, this corresponds roughly to the 
range r = 1.05 fm [20| of the finite-range term in the 
BCP functional. 

For the sake of simplicity, the spin-orbit and Coulomb 
contributions, E s -°- and Eq, will be neglected in the 
present work. Note that, according to the Kohn-Sham 
energy density functional theory, the BCP energy func- 
tional docs not introduce an effective mass m* in the 
kinetic energy part. Hence, the mass m which appears 
in the hamiltonian h is the free nucleon mass. 

With the present prescription to calculate the force on 
a test particle, it is straight-forward to show that the 
total energy defined by 



E tl 



1 MA r? 

_ Ei 



+ E£b\Pp, Pn] 



(13) 



is exactly conserved during the time evolution. 



D. Ground state initialization 

For a given initial distribution function /(r,p,i = to), 
the time evolution at t > to is completely determined 
by Eq. ([3]). Here, we choose as initial state an isolated 
nucleus at rest in its ground state. An obvious require- 
ment for the ground state is that it must be stationary. 
In the present framework, this means that the ground 
state must be calculated within the TF approximation, 
/(r,p) = 0(p — h{r, p)), which is stationary under the 
Vlasov equation (|3]). Here, 9 is the step function and p is 
the chemical potential (Fermi energy). For consistency, 
since the numerical simulation of the Vlasov equation 
requires to smooth the densities and mean fields with 
Gaussians, Eqs. (O and (|9]), the same Gaussians should 
be included in the calculation of the TF ground state 

M- 

Written explicitly, the equations which have to be 
solved sclf-consistcntly are 



P(r) = 
P(r) 



[2m(n - L>(r))] 3 / 2 



3ir 2 h 3 
d 3 s 
WW 



9(jm - U(r)) 
/d2 p(r-s), 



(14) 
(15) 



where U is computed from p according to Eqs. (|12[) and 
(|9]). The chemical potentials (Fermi energies) p p and p n 
are determined from the conditions 



Z = J d 3 r Pp (r), N = J d 3 rp n 



(r) 



(16) 
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FIG. 1: Ground-state densies of neutrons and protons in units 
of po = 0.17 fm" 3 in 22 (a) and 132 Sn (b). Solid lines: self- 
consistent TF densities p(r), dots: corresponding smoothed 
densities p(r) according to Eq. (|15[l . dashes: smoothed densi- 
ties p(r) after a simulation time of t — 2000 fm/c. 



For illustration, the density distributions of protons 
and neutrons in 22 and 132 Sn obtained in this way are 
displayed in Fig. [1] The neutron skin is clearly visible in 
both cases. Note that the TF densities (solid lines) vanish 
at the classical turning points which arc determined by 
fi a = U a . The finite surface thickness stems from the 
self-consistent solution of the TF equations as described 
by Eqs. (O, (15]), (H3), and © with d = 0.7 fm. 

Once the self-consistent TF density distributions are 
obtained, the test particle positions rt are initialized ran- 
domly according to a probability density P(r^) oc p(ri). 
Then, the momenta p i are initialized randomly in a 
sphere with radius PF^i) = ^(37r 2 p(r i )) 1 / 3 in order to 
correctly describe the Fermi motion. Remember that, if 
the Pauli principle is satisfied initially, it is preserved by 
the Vlasov dynamics due to Liouville's theorem [22| . 



E. Numerical parameters and stability 

As mentioned before, the width d of the Gaussians has 
two effects: first, it is necessary to obtain a well-defined 
density distribution, and second, it induces effectively a 
finite-range in the interaction. Most of the results to be 
presented in this paper were obtained with d = 0.7 fm, 
leading to a reasonable smoothing of the mean field as 
discussed in Sees. Ill CI and III Dt This defines the mini- 
mum number of test particles to be used, as there must 
be a sufficiently large number of test particles per vol- 
ume d 3 , otherwise the statistical fluctuations become too 
strong. Here, Af = 2000 test particles per particle were 
used. The equations of motion were solved with the ve- 
locity Verlet algorithm [25| using a time step of 0.1 fm/c. 
After each time step, the mean field U was updated and 
stored in a three dimensional grid with spacing 0.4 fm. 

With these parameters, it was possible to ensure the 
stability of nuclei with large neutron excess, i.e., with 
very weakly bound neutrons in the surface, during the 
long simulation time of 2000 fm/c which is necessary for 
the calculation of the response function (see next subsec- 
tion). The numerical losses due to test particles which 



escape from the nucleus are < 1.3 neutrons in the case of 
22 and < 4 neutrons in the case of 132 Sn. In order to 
illustrate the stationarity of the ground state, we display 
in Fig. [T] the angle-averaged (see appendix) smoothed 
densities p corresponding to the test-particle distribution 
after the simulation (dashes) , which agree very well with 
those calculated during the initialization (dots). During 
the simulation, the kinetic energy (including that car- 
ried away by the test particles escaping from the nucleus) 
drops by ~ 5% in the case of 22 O and by ~ 2% in the 
case of 132 Sn. The relative variation of the total (kinetic 
plus interaction) energy is of the order of 10~ 6 . 

In order to check that the results do not depend sen- 
sitively on the value of the width parameter d, some cal- 
culations with d = 0.5 fm were performed. In this case, 
the number of test particles per particle was increased to 
Af = 5000 in order to limit statistical fluctuations, and 
the spacing of the grid for the mean field was reduced to 
0.3 fm. 



F. Calculation of the response function 

Within the present approach, the collective modes are 
described as time-dependent oscillations of the nucleus 
after a perturbation of the ground state. In order to 
relate these oscillations to the usual response function, 
which is defined in terms of transition probabilities from 
the ground state to excited states, let us temporarily 
leave the semiclassical framework and return to quan- 
tum mechanics. We consider a perturbation hamiltonian 
of the form H eyi {t) = XQS(t), where Q is the excitation 
operator we want to study and A is supposed to be small. 
Then, within linear response theory [27J , the expectation 
value of the operator Q as a function of time is given by 

S(Q)(t) = (Q)(t)-(0\Q\0) 

= -^EK/M°)l 2 -^S^' (17) 
/ 

where |0) is the ground state (of the unperturbed hamil- 
tonian H), |/) is an excited state, and Eq and Ef are the 
corresponding energies. Defining the strength function 
as usual by 



S(E)=Y / \(f\Q\0)\ 2 6(E- 
f 



Ef+Eo), (18) 



we can obtain it from 5 (Q) (t) via a Fourier transform 

1 f°° Et 

S(E) = dtS(Q) (t) sin — . (19) 
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Let us now return to the semiclassical framework. Un- 
der the assumption that Q is a one-body operator, i.e. 



A 



(20) 
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one can calculate its expectation value as 

(Q) (*) = I d\ d 3 p /(r, p, t)q(r, p) , (21) 



where q(r, p) is the Wigner transform of q, which can 
be obtained (at least to leading order in the ft expansion 
[HJ) by replacing the operators f and p in q by their 
classical counterparts r and p. In terms of the test par- 
ticle positions and momenta, this expectation value can 
be expressed as 



AfA 



(22) 



The last point which remains to be explained is how 
the delta function perturbation at t = changes the ini- 
tial distribution function, i.e., in our simulation, the dis- 
tribution of test particles. In principle, one has to solve 
the classical equations of motion J5]) with the perturbed 
hamiltonian ft + Xq5(t) instead of ft. Replacing the delta 
function by a short pulse of length St, one can show that 
in the limit St — > and to leading order in A, the effect of 
the perturbation is to change the positions and momenta 
of the test particles as follows 2 : 

r^ + A^p), p^ Pl -A%^. (23) 
dpi dr, 

To summarize the procedure: First, the test-particle 
distribution at t = is initialized as explained in Sec. 
Ill Dl Then the test-particle positions and momenta are 
changed according to Eq. (|23]l and the mean field U is 
recalculated if necessary (if the excitation operator q de- 
pends on p). After that, the equations of motion © 
are solved simultaneously for all test particles, and the 
mean field U is updated after each time step. In this way, 
one obtains the expectation value (Q) (t) as a function of 
time, and its Fourier transform (|19p gives the strength 
function S(E). 

In practice, it is of course impossible to run the sim- 
ulation to t — oo. Here, the simulations will be stopped 
at imax = 2000 fm/c. In order to avoid oscillations in the 
Fourier transform related to the cut at t max , the strength 
function S(E) will be folded with a Lorentzian of width 
7 = 0.5 MeV, which is equivalent to multiplying the sine 
function in Eq. (19]) by e -T*/ 2ft . 



G. Transition densities and velocity fields 

The delta function perturbation at t = excites simul- 
taneously all modes which can be excited by the operator 
Q. It is therefore difficult to extract the transition density 



2 If q depends only on r or only on p, this result is valid to all 
orders in A. 



and velocity field corresponding to one particular mode. 
What can be done is to calculate the (smoothed) density 
distributions p(r, t) and velocities 



v(r,t) 



1 



MA 

E 



e (r-r t (t)r/d 2 



p(r,t) p(r,t)£< m M(yftdf 



(24) 



(see appendix for more details) as functions of time. 

In the case of a time-even excitation operator, i.e., 
q(r, p) = q(r, — p), the particles get a kick in momentum 
space and the oscillation starts with maximum velocity, 
while the density is not changed at t = 0. The situation 
is opposite if the excitation operator is time-odd, i.e., 
<z(r, p) = —q(r, — p): In this case, the velocity is zero at 
t — 0, while the density is immediately changed due to 
the displacement of the particles in coordinate space. 

In order to find the contribution of a given mode to the 
density and velocity oscillations, one has to choose the en- 
ergy E corresponding to a peak in the strength function 
and compute the transition densities and velocity fields 
as a Fourier transform of p(r, t) and v(r, t), respectively. 
In the case of a time-even excitation operator, one has to 
use 

Sp(r,E)(x / dtp(r,t)sm — , (25) 
Jo h 

f°° Et 
v(r,E)oc / dt-v (r,t) cos — , (26) 
Jo " 

whereas in the case of a time-odd operator, the sine and 
cosine functions must be interchanged. In practice, since 
the simulation runs only to t max = 2000 fm/c, the sine 
and cosine functions are multiplied by an exponential 



damping factor e 



-~it/2h 



7 = 0.5 MeV, as in the strength 



function S{E). 

It should be noted that even if the energy E corre- 
sponds to a peak in S(E), the transition densities and 
velocity fields obtained with this method may still con- 
tain contributions from other modes if those have a width 
which makes their spectrum extend to energy E. 



III. RESULTS FOR THE PYGMY RESONANCE 

A. Electric dipole response 

Since the pygmy resonance is often studied in (7,7') 
experiments, we consider as excitation operator the elec- 
tric dipole operator [28| 



for protons, 
for neutrons, 



(27) 



which is defined such that the center of mass of the nu- 
cleus stays at rest. For the parameter A multiplying the 
operator q, the value A = 25 MeV/c is chosen as a com- 
promise to excite an oscillation which is much larger than 



6 



w 



- 

A 

T ' i 

- ' 1 V\ 

1 .1 ,w \ 

S '<' << 
r 1 : ;i ■ 


i nn 

Sn 

,,,Sn 

132 Sn . 

Lv y../ 'A Iaaa^ : v^~a~v^ 


w I' / 






'A' V / 1 * ;V ' 














r 









200 400 600 



t (fm/c) 

FIG. 2: Electric dipole moment of 100 Sn, lw Sn, and 132 Sn af- 
ter a perturbation with H ayi (t) = AQ<5(i), Q being the electric 
dipole operator and A = 25 MeV/c. 
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FIG. 3: Electric dipole strength of 100 Sn, lie Sn, and 132 Sn. 

the numerical noise due to the finite number of test par- 
ticles, but still small enough so that nonlinearities do not 
play a role. Calculations for different nuclei from oxygen 
isotopes up to 208 Pb were performed. The existence of 
the PDR as a small enhancement in the strength func- 
tion well below the energy of the GDR turned out to be 
a general property of N > Z nuclei, while it is absent in 
N = Z nuclei. 

As a first example, let us discuss the results obtained 
for the three tin isotopes 100 Sn, 116 Sn, and 132 Sn. In 
all three cases, after initializing and exciting the nucleus, 
one observes a damped oscillation of (Q) (t) with the fre- 
quency of the GDR, see Fig. [21 The curves look quali- 
tatively similar for all three nuclei, and in order to see 
anything else than the GDR, one has to look at their 
Fourier transforms. In Fig. [3[ the corresponding electric 
dipole strengths are displayed. In all three nuclei, there 
is a very strong peak at around 14.5 MeV ( 132 Sn) to 16.7 
MeV ( 100 Sn) corresponding to the GDR. The energies of 
the GDR are in reasonable agreement with experimental 
values, whereas the widths are much too small. For in- 
stance, the curve shown for 116 Sn is very well fitted by 



a Lorentzian with energy 15.8 MeV and width 2.4 MeV 
(which includes the artificial width 7 = 0.5 MeV men- 
tioned in the end of Sec. Ill Fft . whereas the corresponding 
experimental energy and width are 15.67 and 4.19 MeV, 
respectively [29| . 

It is in fact not surprising that the widths are too small. 
Since in the linear regime the Vlasov equation is a semi- 
classical version of the RPA fl6l . , the only damping 
mechanism that is present here is Landau damping. The 
fact that in an anharmonic potential different classical 
orbits have different periodicities leads to effects which 
are completely analogous to the splitting of the quantum 
mechanical single-particle levels, as already noticed in 
Refs. (l6l. [l7T|. The main difference to quantum mechani- 
cal RPA calculations is that within Vlasov dynamics the 
strength is not fragmented into many discrete states, but 
it is continuous 3 . The fragmentation due to the cou- 
pling to more complex states like two-phonon or two- 
particle-two-hole states is missing here, as it is in RPA 
[30]. In the semiclassical framework, effects analogous 
to two-particle-two-hole excitations can be included via 
a collision term [TtJ , but this is beyond the scope of the 
present work. 

Let us return to the discussion of the results shown in 
Fig. [3] In the case of 116 Sn, one can see a small amount 
of dipole strength below the GDR which is absent in the 
N = Z nucleus 100 Sn and which becomes a well defined 
peak at 8.6 MeV in the case of 132 Sn. This peak corre- 
sponds to the PDR. For comparison, in experiment, it 
was seen in 130 Sn and 132 Sn at a slightly higher energy of 
approximately 9.8 MeV [3l|. Subtracting the tail of the 
GDR (assuming that it has the same shape as in 100 Sn), 
one finds that in 132 Sn the PDR contributes about 4% to 
the energy- weighted sum rule (EWSR), which happens 
to be in perfect agreement with the experimental value 
from Ref. [3l| . However, the experimental number for 
130 Sn is larger (7% of the EWSR) than that for 132 Sn al- 
though the neutron excess is smaller, certainly due to the 
doubly magic nature of 132 Sn. In a theory without shell 
effects, the results for 130 Sn and 132 Sn are of course al- 
most identical. Nevertheless this comparison shows that, 
in spite of its crudeness, the semiclassical approach is 
capable of giving the right order of magnitude for the 
transition strength. 

In order to see how sensitive these results are to the 
choice of the width parameter d, the calculation for 132 Sn 
was repeated with d = 0.5 instead of 0.7 fm. The main 
effect of this change on the electric dipole response is that 
the GDR is slightly shifted from 14.5 to 15.3 MeV. The 
position of the PDR is almost not affected (the maximum 
of the peak is shifted from 8.6 to 8.7 MeV). The height of 



3 Note that in Refs. [l6l , ll7l the angular momentum of the classical 
orbits was artificially quantized in order to simplify the practical 
calculations and to obtain a discrete spectrum as in RPA. In 
the present work, the angular momentum of the test particles is 
arbitrary, which results in a continuous spectrum. 
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FIG. 4: Electric dipole strength of 16 0, 18 0, 20 O, and 



the peak corresponding to the PDR is slightly reduced, 
but subtracting the tail of the GDR, which is now further 
apart, one finds again that it contributes about 4% of 
the EWSR. The shape of the transition densities and 
velocity fields (see next subsection) of the two modes are 
not changed either, except that they go to zero more 
rapidly at the surface. The conclusion is that the results 
do not depend strongly on the parameter d, and from now 
on we will keep d = 0.7 fm which is the value motivated 
in Sec. HTCl 

As a second example let us consider the even oxygen 
isotopes from 16 O to 22 O. The corresponding strength 
functions are displayed in Fig. As in the case of tin 
isotopes, it can be seen that with increasing neutron ex- 
cess some strength builds up at low energies, which is 
clearly separated from the GDR. Quantitatively, if inte- 
grated up to 15 MeV, its contribution to the EWSR is 
0% for 18 0, 4% for 20 O, and 8% for 22 0. This has to 
be compared with the corresponding experimental num- 
bers which arc 8% for 18 0, 12% for 20 O, and 7% for 22 
[32l ]. It is interesting to notice that, as in the case of 
130 Sn and 132 Sn, the experimental results for the con- 
tribution of the low-lying strength to the EWSR do not 
increase with increasing neutron excess. This must be 
related to shell effects and cannot be reproduced within 
the semiclassical framework. Unlike in 132 Sn, the low- 
lying strength distribution in 22 O is completely spread 
and does not have a clear peak. In this case, it does not 
seem to be appropriate to speak of the pygmy resonance 
as a collective mode. 

Finally, let us discuss some numbers for the nucleus 
208 Pb. The response (not shown) looks qualitatively sim- 
ilar to the 132 Sn case: the pygmy resonance shows up as 
a well-defined peak. This peak is situated at 7.6 MeV, 
while the experimental spectrum has two groups of tran- 
sitions around 5.3 and 7.3 MeV Q. The result for the 
total strength B(E1; f) integrated up to 8 MeV is 2 e 2 fm 2 
within Vlasov and 1.32 e 2 fm 2 in the experiment Q. 
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B. Velocity fields and transition densities 

In order to study the nature of the collective modes, it 
is useful to look at the transition densities and the veloc- 
ity fields. The graphical representation of the transition 
density can be simplified by assuming that the amplitude 
of the oscillation is weak (linear response regime). In this 
case, the spherical symmetry of the ground state and 
the dipole form of the excitation operator imply that the 
transition density can be written as Sp(r) = 8p(r) cos 9, 
where r = |r| and cos 8 = z/r. 

In order to test the calculation of transition densities 
and velocity fields, let us start with a simple example, 
namely with the GDR in the symmetric nucleus 100 Sn. 
The results are shown in Fig. [51 As expected, protons 
and neutrons move against each other in z direction. The 
velocity (Fig. [Sb) is not constant but decreases with in- 
creasing r and gets curved, almost as in the Steinwedel- 
Jensen model in which the radial component of the veloc- 
ity field vanishes at the surface [l4|. Since the Coulomb 
interaction is not included in the present calculation, the 
transition densities and velocity fields of neutrons and 
protons should be exactly opposite to each other. For 
the velocity fields (Fig. [5Jd), this seems to be the case, 
but in the transition densities (Fig. [5jt) a discrepancy is 
present at small radii (< 2 fm). This is clearly a numer- 
ical error. The reason is that at small radii, the angle 
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FIG. 6: Same as Fig. but for the GDR in 132 Sn. 



FIG. 7: Same as Fig. but for the PDR in 132 Sn. 



averaging, which is implicit in the computation of the 
radial functions 5p{r), is less effective in reducing statis- 
tical fluctuations than at large radii. Since the amplitude 
of the oscillation is very small, already small statistical 
fluctuations of the density can lead to an erroneous re- 
sult for the transition density. This is why the transition 
densities at r < 2 fm cannot be trusted. 

After this word of caution, let us look at the more in- 
teresting case of the neutron rich nucleus 132 Sn. Since 
in this nucleus the neutron and proton density distri- 
butions in the ground state are different, one does not 
expect any more that the transition densities and veloc- 
ity fields of neutrons and protons are exactly opposite to 
each other. Generally speaking, in TV =^ Z nuclei, even in 
the absence of Coulomb interaction, the collective modes 
are not exactly isovector or isoscalar ones, but they have 
both isovector and isoscalar components. Let us first 
discuss the GDR which is displayed in Fig. [6] Since the 
transition densities for r < 2 fm are not reliable, the node 
of 5p p (dashed line in Fig. [6^,) is most likely a numerical 
error. Beyond that radius, the shape of the transition 
densities is typical for the GDR. As a consequence of the 
neutron skin, the transition density of neutrons extends 
to larger radii than that of protons. The velocity fields 
(Fig. [Bp) are more surprising: While the proton velocity 
is very similar to the one in 100 Sn, the neutron velocity 
is very different. It seems that in 132 Sn the neutron ve- 
locity is strongly suppressed in the center and enhanced 
in the neutron skin. The origin of this phenomenon is 



not completely understood, but a possible explanation 
could be the coupling between the isovector GDR and 
the isoscalar torus mode, see next section. 

In Fig. [3 the transition densities and velocity fields 
corresponding to the pygmy mode are displayed. As one 
can see, protons and neutrons oscillate mainly in phase. 
The isovector component of the pygmy mode comes from 
the different transition densities in the region of the neu- 
tron skin. The velocity field has a toroidal shape, very 
different from the giant resonance. Such a shape was al- 
ready found in quantum mechanical calculations of the 
velocity field of the pygmy mode in 122 Zr Q and in 208 Pb 
Q . In the literature, it is often said that the PDR is an 
oscillation of the neutron skin against the N = Z core 
of the nucleus [l[ . Due to the toroidal form of the veloc- 
ity field, the neutrons in the neutron skin indeed move 
against the neutrons in the core. However, the image of 
the skin oscillating as a whole against an inert core seems 
to be over-simplified, since the protons have a toroidal 
flow-pattern, too. 



IV. RESULTS FOR THE TORUS MODE 

A. Toroidal excitation spectrum 

Motivated by the toroidal shape of the velocity field of 
the PDR, let us have a closer look at the isoscalar torus 
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toroidal dipole operator, Eq 



mode. Its excitation operator is given by 

?(r,p) = (2r 2 -l(r 2 )) Pz -(p-r)z. 



(28) 



Here, the term cx (r 2 )p z has been added to the operator 
given in Ref. [ll[ in order to make sure that the excitation 
does not displace the center of mass of the nucleus. Note 
that, according to Eq. (|2"3")l . this excitation operator leads 
to a change of both positions and momenta of the test 
particles at t = since it depends on r and p. 

The corresponding strength functions of 100 Sn and 
132 Sn are shown in Fig. [SJ We see that the strength 
is split into two regions below and above ~ 20 MeV. 
The first region contains two peaks at 10 and 13 MeV 
in the case of 100 Sn and three peaks at 8.9, 11, and 14.6 
MeV in the case of 132 Sn. In the second region, there is 
an isolated peak at 32.8 MeV in the case of 100 Sn and 
28.9 MeV in the case of 132 Sn. The nature of the dif- 
ferent modes will be clarified by the analysis of the cor- 
responding transition densities and velocity fields. It is 
interesting to notice that the positions of the modes in 
100 Sn at 10 and 32.8 MeV are in good agreement with 
recent RPA results obtained with the UCOM interaction 
[H[. The mode at 13 MeV corresponds probably to the 
fragmented strength concentrated around 15 MeV in the 
RPA response, see Fig. 1 of Ref. [33j | . 



B. Velocity fields and transition densities 

In Figs. [9lfT2l the transition densities and velocity fields 
corresponding to the peaks of the toroidal dipole response 
of 132 Sn at 8.9, 11, 14.6, and 28.9 MeV arc shown. Let 
us first compare the results for the modes at 8.9 (Fig. 
U) and 11 MeV (Fig. [10]). At first glance, the velocity 
fields look similar for these two modes, but the transi- 
tion densities are completely different. Both modes are 
essentially isoscalar. Comparing the results for the 8.9 
MeV mode in Fig. U with the results for the PDR at 8.6 
MeV in Fig. [7J we see a striking similarity. One can say 
that these two modes are in fact one and the same, only 
excited in two different ways. Note that the lines with 
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FIG. 9: Same as Fig. [5l but for the mode excited by the 
toroidal dipole operator at 8.9 MeV in 132 Sn. 
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FIG. 10: Same as Fig. [5] but for the mode excited by the 
toroidal dipole operator at 11 MeV in 132 Sn. 
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FIG. 13: Same as Fig.0 but for the toroidal mode at 10 MeV 
in 100 Sn. 
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FIG. 12: Same as Fig. [5] but for the compressional dipole 
mode at 28.9 MeV in 132 Sn. 



zero velocity, around which the protons and neutrons cir- 
culate, are at ~ 3 fm and ~ 4 fm from the center of the 
nucleus, respectively, i.e., they lie inside the core of the 
nucleus and not in the surface. The velocity field of the 
mode at 11 MeV looks even more like a torus, since the 
shape of the velocity field is more rounded than in the 
mode at 8.9 MeV. What about the third mode at 14.6 
MeV (Fig.rTTj)? First of all, from the transition densities, 
one sees that this mode has mainly isovector character. 
Comparing with the results for the GDR at 14.5 MeV in 
Fig- El one concludes that the mode at 14.6 MeV is in fact 
the GDR, which is excited by the isoscalar toroidal dipole 
operator due to the strong neutron excess in 132 Sn. Fi- 
nally, the high-lying isoscalar dipole mode at 28.9 MeV 
(Fig. [T2"|) has a completely different nature. As can be 
seen from the velocity field, this mainly isoscalar mode 
exhibits a compressional motion, and for this reason it is 
usually called the compressional dipole mode. 

In order to get a better understanding of the two low- 
lying modes at 8.9 and 11 MeV, let us look at the corre- 
sponding modes of 100 Sn which lie at 10 and 13 MeV, sec 
Figs. [13] and [H Wc observe that the mode at 10 MeV 
in 100 Sn (Fig. HI?]) has qualitatively the same velocity field 
and transition density as the mode corresponding to the 
PDR in 132 Sn (Fig. [9]). From this one may conclude that 
the existence of this mode does not require the presence 
of a neutron skin. The neutron excess is only needed 
in order to be able to probe this mode with the electric 
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FIG. 14: Same as Fig.0 but for the toroidal mode at 13 MeV 
in 100 Sn. 



dipolc operator. Looking closely at the velocity field of 
this mode (Figs. l9l and fT3|). one realizes that the veloc- 
ity field is almost constant in the center of the nucleus, 
contrary to the velocity field of the toroidal mode which 
lies at slightly higher energy [11 MeV in 132 Sn (Fig. [I"0|) 
and 13 MeV in 100 Sn (Fig. Q3J}, respectively]. This can 
be interpreted in the sense hat the lower one of the two 
modes is an oscillation of the surface (not necessarily the 
neutron skin, since 100 Sn does not have one) against the 
core, whereas the higher one is the original torus mode 
which, in the framework of nuclear fluid dynamics, exists 
already in a uniform sphere. This interpretation is cor- 
roborated by the transition densities, which in the case 
of the higher mode (Figs. [TUl and [T4|) are much more con- 
centrated in the inner part of the nucleus, while those of 
the lower mode (Figs. l9l and IT3|) are much stronger in the 
surface region. Another support for this interpretation is 
the energy of this mode: According to Ref. [l2j , in a uni- 
form sphere and within nuclear fluid dynamics, the torus 
mode should lie at 65 - 85A" 1 / 3 MeV, i.e., at 14 - 18.3 
MeV in the case of 100 Sn and 12.8-16.7 MeV in the case 
of 132 Sn. This is slightly higher, but not very far from 
the modes found here at 13 and 11 MeV, respectively. 



V. SUMMARY AND DISCUSSION 

In this paper, electric and isoscalar dipolc excitations 
were studied within the semiclassical TF plus Vlasov ap- 
proach. As interaction, the bulk part of the BCP func- 
tional was employed, which was smoothed out in space 
in order to mimic the effect of the neglected finite-range 
term. The TF equation for the ground state as well as 
the Vlasov equation for the dynamics were solved numeri- 
cally without any further simplifying assumptions. Com- 
pared to fully quantum mechanical Hartrcc-Fock plus 
RPA calculations, the present method is missing the shell 
effects, what can be useful if one is interested in generic 
properties and average trends. 

The electric dipole response of oxygen and tin isotopes 
was discussed (calculations for other nuclei like calcium 
and lead were performed but not shown). In all cases, 
low-lying strength corresponding to the PDR was found 
in the very neutron-rich isotopes, however in the neutron- 
rich oxygen isotopes the strength was spread over a large 
energy range, in accordance with the strong fragmenta- 
tion of the strength in quantum mechanical calculations 
@ , so that one cannot speak of a collective mode in this 
case. This shows that the existence of the PDR is a 
generic property of neutron-rich nuclei and does not rely 
on a specific structure of single-particle levels. In heavier 
nuclei, the PDR is found to be a collective excitation. 
The obtained energies and transition probabilities are 
roughly in agreement with experimental data (at least 
as well as it can be expected in a theory without shell ef- 
fects). The transition densities and velocity fields of the 
pygmy mode were analysed and it was found that the 
velocity field has a toroidal shape, in agreement with the 
findings of earlier quantum mechanical calculations [H, H| . 
The vortex line lies inside the core of the nucleus, and the 
torus can therefore be seen in both neutron and proton 
velocity fields, which suggests that the popular picture of 
the neutron skin oscillating against a static N = Z core 
is oversimplified. 

In order to compare the PDR and the torus mode in 
detail, the response to the isoscalar toroidal operator was 
studied. In the example of 132 Sn it was found that, due to 
the mixing of isoscalar and isovector modes, the predom- 
inantly isoscalar PDR and the predominantly isovector 
GDR can be seen in both the El and the toroidal strength 
functions. In addition, the toroidal response exhibits two 
peaks which do not show up in the electric dipole re- 
sponse: a second toroidal mode which lies slightly above 
the PDR, and a compressional dipole mode. The two 
toroidal modes exist also in the N = Z nucleus 100 Sn. 
Hence, the existence of these modes, including the lower 
one which in the case of 132 Sn was identified with the 
PDR, does not rely on the presence of a neutron skin. 
If this interpretation is correct, the reason why the PDR 
is only seen in nuclei with large neutron excess is simply 
that otherwise the El strength of this mode is too small 
to be seen. However, it was recently pointed out that, 
if the isospin symmetry breaking effect of the Coulomb 
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interaction is taken into account, a small contribution of 
this mode can be seen even in the electric dipole response 
of N = Z nuclei [33|). Another conclusion which can be 
drawn from this result is that the PDR, like the torus 
mode, cannot be described in a hydrodynamical picture, 
but its existence relies on the "elasticity" of the nuclear 
medium due to Fermi-surface deformation. This fact was 
also stressed in Rcf. 
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The two toroidal modes are apparently qualitatively 
different, although their velocity fields look quite simi- 
lar: The higher mode corresponds to the torus mode of 
an elastic sphere, which has been discussed in the litera- 
ture for many years [lol [l2j , whereas the lower one cor- 
responds more to an oscillation of the core against the 
surface (but not necessarily against the neutron skin) , 
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qualitatively similar to the modes discussed in Ref. 

Of course, the present approach has some shortcomings 
and is not meant to replace more sophisticated quantum- 
mechanical calculations. Since it is a semiclassical for- 
malism, shell effects cannot be described. On the one 
hand, this results in clear pictures for the different modes, 
but on the other hand, it is of course a simplification 
which makes the detailed comparison with experiment 
and with more realistic calculations difficult. In addi- 
tion, a couple of approximations were made which give 
probably rise to systematic deviations. For example, the 
Coulomb interaction was omitted. One would also ex- 
pect an important effect from pairing, since it affects 
in particular rotational motion and excitations involving 
Fermi surface deformation. Nevertheless, the obtained 
results are surprisingly reasonable and will maybe serve 
as a motivation for a more detailed study of the relation- 
ship between the pygmy and the torus mode within fully 
quantum mechanical approaches. 



The dipole excitations considered in this work destroy 
the spherical symmetry, but not the cylindrical symmetry 
with respect to the z axis. In Figs. [5][7] and IMT4"1 this 
symmetry was used to reduce the statistical fluctuations 
by averaging the density p and the components j±_ and j z 
of the current density over the azimuthal angle cj> (j^ = 
for the excitation operators under consideration). After 
this averaging, the final expressions for p and j in terms 
of the test particle positions and momenta read: 



MA _ (— ' O'+fx+^Xi 



, (A2) 



MA - ( " i) y +rL 2r r 



MA 



j±(r±,z) = ^P±i cos0 ri 



2,2,2 
+ ' ±+' ±j 

3 flrxr 
h 



i=l 



N{^df 



\ d 2 J 
(A4) 



where Iq and I\ are modified Bessel functions [35| and 
4>npi is the difference of the azimuthal angles of and 
Pz, i-e-, 



P±i COS 4> nPi = 



XiPx 



UiPyi 



V x i+y'i 



(A5) 



Appendix B: Transition probabilities and energy 
weighted sum rule 
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Appendix A: Angle averaged densities and velocity 
fields 



All simulations were done in three dimensions with- 
out any imposed symmetries. However, for the graphical 
representation of the results it is advantageous to average 
the densities and velocity fields over the angle in order 
to reduce the statistical noise due to the finite number of 
test particles. 

Let us start with the ground state density distribu- 
tions shown in Fig. [TJ Within the TF approximation, 
the ground states of all nuclei are spherical. Therefore, 
the ground state densities p were averaged over the full 
solid angle. In terms of the test-particle positions r^, the 
angle-averaged densities can be expressed as 
MA e _ (r _ ri) 2 /d 2 _ e _ (r+ri) 2 /d 2 

fa) = E M^rrA ■ (A1) 



Often, the B(E\) value of the pygmy mode or its con- 
tribution to the EWSR are used as a measure for the 
strength of the py gmy mode. According to the defini- 
tions given in Ref. [28[, the reduced transition probabil- 
ity B(E1;0 — > 1) from the I\ = ground state to a 
I2 = 1 excited state can be related to the strength func- 
tion S(E) corresponding to the electric dipole operator 
([2TJ) as follows: 



dB(El;0^1) 9e 2 
dE = 4^ 5(E) 



(Bl) 



Since the BCP functional does not introduce an effective 
mass (m* — to), the EWSR (Thomas- Rciche-Kuhn sum 
rule) 



/>oo 

/ dEES(E) = 
Jo 



h 2 NZ 
2m A 



(B2) 



should be exactly fulfilled. In the numerical results dis- 
cussed in Sec. IIII A[ the deviation from the exact result 
is less than 1%. 
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